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Abstract 

The most efficient MC weights for the calculation of physical, canonical expectation values are not necessarily 
those of the canonical ensemble. The use of suitably generalized ensembles can lead to a much faster convergence 
of the simulation. Although not realized by nature, these ensembles can be implemented on computers. In recent 
years generalized ensembles have in particular been studied for the simulation of complex systems. For these 
systems it is typical that conflicting constraints lead to free energy barriers, which fragment the configuration 
space. Examples of major interest are spin glasses and proteins. In my overview I first comment on the strengths 
and weaknesses of a few major approaches, multicanonical simulations, transition variable methods, and parallel 
tempering. Subsequently, two applications are presented: a new analysis of the Parisi overlap distribution for the 
3d Edwards- Anderson Ising spin glass and the helix-coil transition of amino-acid homo-oligomers. 

Key words: Markov chain Monte Carlo simulations, multicanonical, transition matrix, parallel tempering, first order 
phase transitions, complex systems, spin glasses, proteins, peptides. 
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1. Introduction 

Markov chain Monte Carlo (MC) simulations are 
an indispensable tool for the study of physical mod- 
els. To calculate expectation values in the Gibbs 
canonical ensemble, most simulations use Boltz- 
mann weights. However, in the course of time it 
has become clear that the Boltzmann factor is not 
always the most efficient weight to generate the de- 
sired configurations of the canonical ensemble. In 
particular before embarking on a large-scale com- 
puter simulation, one of the questions which ought 
to be addressed is "What are suitable weight fac- 
tors for the problem at hand?" Nowadays simula- 
tions of more general ensembles are employed for 
an increasing number of applications in physics, 



chemistry, structural biology and other areas, for 
reviews see [1-4] . 

For instance, the multicanonical (MUCA) 
method [5-10] calculates in one simulation the 
canonical ensemble at many temperatures. Its 
best-established application is calculations of in- 
terface tensions at first order phase transitions. 
There, canonically rare configurations are simply 
enhanced, the statical case. This is similar to the 
importance sampling idea, which for the calcu- 
lation of canonical expectation values led from 
naive sampling to the Markov process sampling 
with the Boltzmann weights. The additional com- 
plication is now that the new weights are a priori 
unknown. The nowadays most important appli- 
cations of the MUCA method are for complex 
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systems [7,8]. There, the rational of MUCA simu- 
lations is to flatten free energy barriers, such that 
the dynamics of the simulation becomes improved. 
The enhancement of particular configurations can 
no longer be explicitly controlled and the issue of 
optimizing MC algorithms for complex systems 
is far from being well understood, although some 
improvements have been quite spectacular. 

In the next section we introduce the MUCA 
method, using as illustration the 2d ten-state Potts 
model with its strong first-order phase transition. 
Subsection 2.1 discusses the weight recursions of 
Rcf. [1,9] and [10]. Random walk and transition 
matrix methods [11-18] generalize the MUCA ap- 
proach and are summarized in section 3. Distinct 
from these methods are multiple Markov chain 
algorithms [19-21] (parallel tempering). They 
are introduced in section 4, where their recent 
combinations [22] with MUCA methods are also 
sketched. In section 5 we turn to applications for 
complex systems. Instead of being exhaustive (see 
the reviews [1-3]) we focus on two recent studies: 
a new analysis of the Parisi overlap distribution 
for the 3d Edwards- Anderson Ising spin glass [23] 
and the helix-coil transition of the polyalanine 
protein [24,26]. Outlook and conclusions are given 
in section 6. 



2. Multicanonical Simulations 

Multicanonical simulations are best established 
for investigations of first-order phase transitions. 
For them pscudocritical points (3 C (L), where L is 
the lattice size, exist such that the canonical en- 
ergy density P(E) — P(E; L) becomes double- 
peaked and the maxima at -E^iax < ^ ax are of 
equal height P max = P(E^ ax ) = P(^ ax ). In be- 
tween these values a minimum is located at some 
energy E m i n . Configurations at E m i n are exponen- 
tially suppressed like 



P min = P(E min ) = c f LP exp(-/M) 



(1) 



where f s is the interface tension and A is the mini- 
mal area between the phases, A = 2L d ~ 1 for an L d 
lattice, Cf and p are constants. To determine the 



interface tension with Binder's histogram method, 
one has to calculate the quantities 



- with R(D = -f-w 



A(L) 



(L) 



(2) 



and to make a finite size scaling (FSS) extrapola- 
tion of f s (L) for L — > oo. However, for large sys- 
tems a canonical MC simulation will practically 
never visit configurations at energy E = E m - ln and 
estimates of the ratio R(L) will be very inaccu- 
rate. The terminology supercritical slowing down 
was coined to characterize thus an exponential de- 
terioration of simulation results with lattice size. 
The MUCA method overcomes this problem by 
sampling, in an appropriate energy range, with an 
approximation 



w mu (k) = w mu (E {k) ) = e 



b(_E< fc »)_E< fc »+a(B lK ') 



(3) 



to the inverse density of states l/n(E^). Here the 
function b(E) is the inverse microcanonical tem- 
perature at energy E and a(E) is the dimension- 
less, microcanonical free energy. One samples in- 
stead of the canonical energy distribution P(E) a 
new MUCA distribution 



Pmu(E^ — C mu 7l(-E/) W mu (E^) ~ C mu . 



(4) 



The desired canonical distribution is obtained 
by rc-wcighting, which is rigorous, because the 
weights w mu (E) used in the actual simulation are 
exactly known. With the approximate relation (4) 
the average number of configurations sampled does 
no longer depend strongly on the energy and accu- 
rate estimates of the ratio (2) R(L) — P m in/-Pmax 
become possible. 

The MUCA method requires two steps 

(i) Obtain a working estimate of the weights (3). 
Working estimate means that the approxi- 
mation to the inverse density of states has to 
be good enough to ensure movement in the 
desired energy range. 

(ii) Perform a Markov chain MC simulation with 
the weights w mu (k). Canonical expectation 
values are found by re- weighting to the Gibbs 
ensemble and standard jackknife methods al- 
low reliable error estimates. 

For the 2d 10-state Potts model figure 1 re- 
produces a thus obtained MUCA histogram to- 
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Fig. 1. The multicanonical energy histog (E) to- 

gether with its canonically re- weighted energy density P(E) 
for the 2d 10-state Potts model on a 70 X 70 lattice [6]. 



gether with its canonically re-weighted energy 
histogram [6]. The canonical energy density of 
this figure gives the estimate 2f s (L) of the inter- 
face tension on the 70 x 70 lattice. Combining 
the such obtained estimates from several lat- 
tices allows the infinite volume FSS extrapolation 
2/ s = 0.0978(8). At the time of its publication 
this number was in strong disagreement with other 
numerical results. Afterwards the exact value was 
discovered to be 2/ s = 0.094701... [27]. This suc- 
cess helped the acceptance of the method, which 
by now has been applied to a large number of 
models [1]. 

There exist many variants of the MUCA 
method. MUCA refers to calculations of canon- 
ical expectation values for a temperature range 
and re-weighting has to be done in the internal 
energy. Similarly, other physical quantities can be 
considered, e.g. multimagnetical [28] refers to sim- 
ulations which give results for a certain range of 
the magnetic field. A variant for cluster updates 
due to Janke and Kappler [29] is called multi- 
bondic. In Ref.[30] the multi-overlap algorithm was 
introduced, which focuses on the Parisi order pa- 
rameter of spin glasses (see section 5). Combining 
MUCA with multigrid methods has been explored 
by Janke and Sauer [31]. For molecular dynamics, 
Langevin and hybrid MC variants see Hansmann 
et al. [32] and, with emphasise on lattice gauge 
theory, Arnold et al. [33]. 



2.1. Multicanonical Recursion 

For spin systems with first-order phase transi- 
tions the FSS behavior is relatively well known. 
Provided the steps between system sizes are not 
too large, it is then possible to obtain working es- 
timates of the w mu (E) weights by means of a FSS 
extrapolation from the already simulated smaller 
systems [6]. However, this method cannot be ap- 
plied to complex systems and to get a working es- 
timate of the MUCA parameters by means of a re- 
cursion is then at the heart of the method. 

As the overall normalization is irrelevant for up- 
dating purposes, it is sufficient to consider nearest 
neighbor (in energy) ratios of the weights (4) 



R n {E) 



w{E) 
w(E + e) 



(5) 



where e is the smallest available energy step and n 
refers to the n th iteration. Initially, w°(E) — 1 is a 
good starting value, because the system can then 
move freely. The final equation for the recursion, 
as derived in [1] , reads 



R n+1 {E) = R n (E) 



H n (E + e) 
H n (E) 



9o(E) 



(6) 



where the exponent g r o(E) is recursively deter- 
mined by the equations 



9o(E) 
9o(E) = 



g n (E)+g%(E) ' 
H n (E + e) H n (E) 



and 



H n (E + e) + H n (E) 
g n+1 (E)=g n (E)+gZ(E), g°(E) = 0. 
Typically, we want to cover an energy range 



(7) 

(8) 
(9) 



The optimum for a flat energy distribution is given 
by a random walk in the energy. This implies a 
CPU time increase ~ V 2 to keep the number of 



E u 



E„ 



£ max transitions constant. The 



recursion (6) needs an additional ~ U 5 (opti- 
mum) number of attempts to cover the entire 
range. For the two-dimensional Ising model on 
L x L lattices, with E min — —2L 2 (the ground 
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Table 1 

Sweeps per i? maI = — > i? min = -L 2 — > £ ma x = 
transition for Lx L Ising models. Here n refers to the first 
transition during the recursion (6), r rec to the subsequent 
transitions during the recursion and r pro£ j to the transitions 
during the production part. 



L 


Tl 


Tree 


T prod 


10 


1654 (66) 


576 (14) 


545 (14) 


20 


19573 (730) 


3872 (120) 


3855 (110) 


30 


79743 (2600) 


12820 (350) 


11428 (300) 


40 


192274 (6900) 


25901 (820) 


27898 (1700) 


50 


472103 (13000) 


48227 (1600) 


45803 (2300) 


60 


810042 (23000) 


78174 (2300) 


78288 (7200) 


80 


2917722 (56000) 


178456 (8900) 


184236 (9100) 



state) and £ max = (the average energy at (3 = 0), 
table 1 gives the number of sweeps per transition. 
The first transition is tedious to find, additional 
transitions are much faster. 

A new, possibly more efficient, recursion was re- 
cently proposed by F. Wang and Landau [10] (my 
interpretation of the merits of the method differs 
somewhat from the authors). They perform up- 
dates with estimators n(E) of the density of states 



p{Ei — > E 2 ) = min 



n(.Ei) 
n{E 2 y 



(10) 



Each time an energy level is visited, the estimator 
is updated according to 



n{E) - n(E) f 



(11) 



where, initially, n(E) = 1 and / = /o = e 1 . Once 
the desired energy range is covered, the factor / is 
refined, 



fl - \JJ , fn+l = \J~Tn. 



+ 1 



(12) 



until some small value like / = e~ 8 = 1.00000001 
is reached. A fast convergence towards a good es- 
timator of the spectral density is reported. After- 
wards the usual, MUCA production runs may be 
carried out. A comparison of the recursions (6) 
and (11,12) is subtle, because in neither case is the 
optimal schedule or termination point known. 



3. Transition Matrix Methods 

The basic idea of the approach [11-18] is to 
employ general transition matrix elements instead 
of just weight factors. An example is Wang's [17] 
random walk algorithm, which uses the transition 
probabilities [11,13-15] 



p(k — > k') = min ( 1 



N(E + AE,-AE) 
N(E,AE) 



(13) 



where N(E, AE) is the microcanonical average for 
the number of transitions from the configuration k 
with energy E to the configuration k! with energy 
E' = E + AE. A flat histogram is obtained and 
the configuration dependence is reduced to energy 
dependence. The density of states is obtained from 
the equation [13-15] 



n(E) N(E, AE) = 



(14) 



n(E + AE) N(E + AE, -AE) . 

The history and general merits of the approach are 
reviewed in [4]. There, it is claimed that, in some 
situations, equation (14) works more efficiently 
than the Wang-Landau recursion (11,12). How- 
ever, the issue may not be settled. In particular, it 
seems [17] that in the random walk MC the use of 
estimators for N(E, AE) faces more serious prob- 
lems than the use of estimators for the weights 
in MUCA simulations. After defining reasonable 
schedules and termination points, one could com- 
pare the transition times of table 1 with those 
achieved with Wang's random walk approach. 



4. Parallel Tempering 

The developments sketched above have to be dis- 
tinguished from the method of multiple Markov 
chains (parallel tempering), which was introduced 
by Geyer [19] and, independently, by Hukusima 
and Nemoto [20]. The latter work was influenced 
by the simulated tempering [34] method, which in 
turn can be understood as a special case of the 
method of expanded ensembles [35] . Parallel tem- 
pering is particularly well suited for parallel pro- 
cessing on clusters of workstations. 
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Parallel tempering performs n canonical MC 
simulations at different /3-values with Boltzmann 
weight factors 

w B ,i{E {k) ) = e-* EW , i=l,...,n (15) 

with /3i < P2 < ••• < Pn-i < fln and allows ex- 
change of neighboring (3- values 

Pi-i < — > Pi for i = 2,...,n. (16) 

These transitions lead to the energy change 

AE = (-A-i^f } - (17) 

which is accepted or rejected according to the 
Metropolis algorithm. The /^-values have to be 
determined such that a reasonably large accep- 
tance rate is obtained for the (3 exchange (16) and 
Rcf. [20] employs a recursive method. 

Remark: The method works for dynamical but 
not for statical supercritical slowing down, because 
each member of the discrete set of weight factors 
samples still a Boltzmann distribution (e.g. it is 
not a valid method for calculating interfacial ten- 
sions) . This can possibly be overcome by applying 
the parallel tempering idea to Gaussian distribu- 
tions instead of the canonical ensemble (communi- 
cated to the author by T. Neuhaus). 

Ref.[22] explores earlier introduced [36] com- 
binations of parallel tempering, called replica- 
exchange method in [22], with MUCA methods. 
The replica-exchange MUCA algorithm performs 
a short replica-exchange simulation to determine 
the MUCA weight factors and continues with 
the production part of a regular MUCA simula- 
tion, whereas MUCA replica-exchange works with 
replicas of the MUCA ensemble in different energy 
ranges. 



5. Applications to Complex Systems 

The relevance of MUCA methods for complex 
systems was pointed out early [7] . In these systems 



one encounters large free energy barriers due to dis- 
order and frustration. Multicanonical simulations 
try to overcome the barriers through excursions 
into the disordered phase. In addition, one may 
try to identify and enhance relevant rare configu- 
rations. Here, we limit our interest to two recent 
studies of (a) spin glasses and (b) proteins. For re- 
views of applications to complex systems see [1-3] . 

5.1. Spin glasses 

The Parisi overlap parameter q is defined as the 
overlap of two replicas of the system at identical 
temperatures 

N 

g = XX-?- (18) 

i=i 

The free-energy structure of the spin glass system 
is thought to be intimately related to the barriers 
in the Parisi order parameter distributions. The 
multi-overlap variant [30] of the MUCA method en- 
hances the minima of the Parisi overlap parameter 
densities. Using the method in large-scale simula- 
tions led to a variety of new insights, for instance 
about the self-averaging properties of free energy 
barriers [37]. The enhancement is most dramatic 
for the tails of the averaged Parisi order parameter 
distribution. Recently, this allowed us to demon- 
strate a conjectured relation between the overlap 
probability density and extreme order statistics 
over about 80 orders of magnitude [23] , see figure 2. 

5.2. Proteins 

Proteins are linear polymers with the 20 natu- 
rally occurring amino acids as monomers. Chains 
smaller than a few tens of amino acids are called 
peptides. The problem is to predict the folded con- 
formation of proteins and peptides solely from their 
amino acid sequence. For many years the emphasis 
of numerical investigations has been on finding the 
global minimum potential energy and the major 
difficulty encountered is the multiple minima prob- 
lem. Molecular dynamics has been the numerical 
method of first choice, but the fraction of stochas- 
tic investigations shows an increasing trend. 
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Fig. 2. Tails of the rescaled overlap (18) probability density 
in comparison with predictions of extreme order statistics 
[Eq.(ll)], see [23] (T = 1.14 is at the freezing temperature). 



The major advantage of MUCA and related 
methods in the context of proteins is that they 
allow for investigations of the thermodynamics of 
the entire free energy landscape of the protein. 
This was realized by Hansmann and Okamoto [8] 
when they introduced MUCA sampling to the 
problem of protein folding and, slightly later, by 
Hao and Scheraga [38] . Since then numerous appli- 
cations have been performed and the simulations 
have been quite successful for peptides. By now a 
quite extensive literature exists, which is compiled 
in [2,3]. In Ref. [39] the weight recursion (6) is 
adapted to these systems with continuous energy. 

A particularly nice application is the helix for- 
mation of polyalanine by means of a MUCA sim- 
ulation [24,25], depicted in figure 3. Up to N = 30 
amino acids are studied numerically and a phase 
transition develops at T c = (541 ± 8) K in the infi- 
nite volume limit (pseudocritical temperature for 
finite systems are in the range from 427 K (N = 
10) to 518 if (N = 30)). FSS estimates of the 
critical indices a, v and 7 are consistent with a 
first-order transition. Progress is made concerning 
the inclusion of aqueous solutions into the simula- 
tion [26]. 




Fig. 3. A helix configuration from a multicanonical simu- 
lation of polyalanine [24] (courtesy Ulrich Hansmann and 
Yuko Okamoto). 

6. Outlook and Conclusions 



Sampling of broad energy distributions allows 
supercritical slowing down to be overcome. This 
is well established for first-order phase transitions. 
Systems with conflicting constraints remain, de- 
spite some progress, notoriously difficult and for 
them most hope lies in achieving further algorith- 
mic improvements. Reasonably large peptides, at 
the border to proteins, are promising systems for 
a complete coverage. Transition matrix extensions 
of the MUCA method and its combinations with 
parallel tempering appear promising. 

This work was, in part, supported by the U.S. 
Department of Energy under contract DE-FG02- 
97ER41022. 
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